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Abstract 

We have calculated the pair correlation functions of a fluid interacting via the Gay-Berne (n- 
6) pair potentials using the Percus-Yevick integral equation theory and have shown how these 
correlations depend on the value of n which measures the sharpness of the repulsive core of the 
pair potential. These results have been used in the density-functional theory to locate the freezing 
transitions of these fluids. We have used two different versions of the theory known as the second- 
order and the modified weighted density-functional theory and examined the freezing of these 
fluids for 8 < n < 30 and in the reduced temperature range lying between 0.65 and 1.25 into 
the nematic and the smectic A phases. For none of these cases smectic A phase was found to be 
stabilized though in some range of temperature for a given n it appeared as a metastable state. 
We have examined the variation of freezing parameters for the isotropic-nematic transition with 
temperature and n. We have also compared our results with simulation results wherever they are 
available. While we find that the density-functional theory is good to study the freezing transitions 
in such fluids the structural parameters found from the Percus-Yevick theory need to be improved 
particularly at high temperatures and lower values of n. 

PACS numbers: 61.30 Cz, 62.20 Di, 61.30Jf 
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I. INTRODUCTION 



In the case of non-spherical molecules the anisotropic nature of the intermolecular inter- 
actions can give rise to new phases (liquid crystals) Q that are absent when simple spherical 
molecules are considered. Depending upon the shape and the size of molecules and upon 
the external parameters (temperature, pressures, etc.) a system may show a wide variety 
of phenomena and transitions in between the isotropic liquid and the crystalline solid. All 
these phases including that of the isotropic liquid and the crystalline solids are characterized 
by the average positions and orientations of molecules and by the intermolecular spatial and 
orientational correlations. The determination of phase diagram of such a system from the 
intermolecular potential is one of the most challenging problems of the statistical mechanics. 

The molecules of systems which exhibit liquid crystalline phases are generally large and 
have group of atoms with their own local features. In general it is difficult to know the 
true nature of the potential energy of interaction between such molecules. Attempts have, 
however, been made to find the potential energy of interactions between two such molecules 
using different approximations. One such method is to sum the interatomic or site-site 
potentials between atoms or between interaction sites. In another and more convenient 
approach one uses rigid molecules approximation in which it is assumed that the inter- 
molecular potential energy depends only on the position of the centre of mass and on their 
orientations. If, however, our interest is to relate the phases formed and their properties to 
the essential molecular factor responsible for the existence of liquid crystals, it is desirable 
to use a phenomenological description, either as a straightforward model unrelated to any 
particular physical systems or as a basis for describing by means of adjustable parameters 
between two molecules. Most commonly used models are hard-ellipsoids of revolution, hard 
spherocylinders @, cut-sphere, the Kihara core model H and the Gay-Berne || model. All 
these are single site models and refer to rigid molecules of cylindrical symmetry. Even for 
these simple models calculating the complete phase diagram is difficult. 

The Gay-Berne potential, in particular, is proving to be a valuable model with which to in- 
vestigate the behavior of liquid-crystals in recent years using computer simulation techniques 
|5|, £| 0| . In this paper we consider a general Gay-Berne (GB) model with n-6 dependence 
on the shifted and scaled separation, R, between the uniaxial particles. 
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where 



u(ei, ej, f) = 4e(ei, ej, r)(iT" - R~ 6 ) 
(r-<r(ei,ej,f) 
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While unit vectors ei, ej indicate the orientations of symmetry axes of particles i and j, 
the orientation of the vector joining them is denoted by the unit vector r. The dependence 
of the contact distance on the orientations of the particles and the interparticle vector is 



<7(ei,ej,r) = a 



i-x 



(ej.f) 2 + (ej.r) 2 - 2x(e i .r)(e j .r)(e i .e j 
1 - x 2 (e;.ej) 2 



(1.3) 



where o"o is the contact distance for the cross configuration (eVej = eYr = ej.r = 0). The 
parameter x is a function of the ratio x (= which is defined in terms of the contact 
distances when the particles are end-to-end (e) and side-by-side (s), 



X 



xj-l 
xl + l 



;i-4) 



This vanishes for a sphere and tends to the limiting value of unity for an infinitely long 
rod. The orientational dependence of the potential well depth is given by a product of two 
functions, 



e(ei, ej, r) = e e l/ (e i , ej)e' M (ei, ej, r) 



1.5) 



where the scaling parameter e is the well depth for the cross configuration. The first of 
these functions 



e(e i ,e j ) = [l- X 2 (e i .e j ) 2 ]^ (1.6) 

clearly favours the parallel alignment of the particles and so aids liquid crystal formation. 
The second function has a form analogous to a(Si, ej, r), i.e. 



e'(ei,ej,r) 



i-x 



e ; .r) 2 + (ej.r) 2 - 2x / (e i .r)(ej.r)(e i .e j 



1 - X /2 (ei.ej) 2 

where the parameter x' is determined by the ratio of the well depths, k'(= ^), via 



X 



;i.7) 



(1.8) 
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The potential contains four parameters (xq, k', /i, v) which determine the anisotropy in the 
repulsive and attractive forces, in addition to two parameters (00, eo) which scale the distance 
and energy, respectively. The ratio of the end-to-end and side- by-side contact distance, x , 
is related to the anisotropy of the repulsive forces and it also determines the difference in 
the depth of the attractive well between the side-by-side and the cross configurations. The 
parameter k' is the ratio of the well depth for the side-by-side and end-to-end configurations. 
While Xo determines the ability of the system to form an orientationally ordered phase, k' 
determines the tendency of the system to form a smectic phase [ff|.The other two parameters 
fj, and v influence nematic and smectic forming character of the anisotropic attractive forces 
in a more subtle way. 

In almost all of the simulation and theoretical studies to date n has been taken equal to 
12. The value of n defines the nature of the repulsion; the higher the value of n the harder 
is the nature of the repulsion. In Fig.l we plot u*(r, Qi, ^2)(= u(r, Q±, Q 2 )/^o) as a function 
of separation for some fixed orientations with n=10 and 18. It shows that as n increases 
the importance of attractive interaction increases for all orientations. In the present paper 
we investigate the effect of variation of n i.e. variation of the range of repulsion on the 
properties of molecular liquids and on its freezing transition. 

The paper is organized as follows: In Sec. II, we describe the solution of the Ornstein- 
Zernike equation using the Percus Yevick closure relation for pair correlation functions. 
Section III discusses the essential details of density functional formalism applied to study 
the freezing of molecular fluids into ordered phases. The results are given and discussed in 
section IV. 

II. PAIR CORRELATION FUNCTIONS: SOLUTION OF THE PERCUS-YEVICK 
EQUATION 

The single particle density distribution p(l) defined as 

P (i) = P (r, n) = (£ 5(v - P ,)*(n - n,)) (2.1) 

1=1 

where and fij give the position and the orientation of i th molecule, the angular bracket 
represents the ensemble average and the 5 the Dirac delta function, is constant independent 
of position and orientation for an isotropic fluid. It therefore contains no information about 
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the structure of the system. The structural information of an isotropic fluid is contained in 
the two-particle density distribution p(l,2) which gives the probability of finding simulta- 
neously a molecule in a volume element dridfli centered at (r 1; fli) and a second molecule 
in a volume element dr 2 dfl 2 centered at (r 2 , f2 2 ). p(l, 2) is defined as 



p(i, 2) = P (n, n i; r 2 , n 2 ) = (£ <j(n - ri )a(«i - n^fo - rj >(ft 2 - n,)) (2.2) 



The pair correlation function g(l, 2) is related to p(l, 2) by the relation 

(2 - 3) 

Since for the isotropic fluid p(l) = p(2) = p/ = where < N > is the average number 
of molecules in volume V, 

^(r,ni,n 2 ) = P (T,n 1 ,n 2 ) (2.4) 

where r = (r 2 — ri). In the isotropic phase p(l, 2) depends only on the distance |r 2 — ri| = 
r, the orientation of molecules with respect to each other and on the direction of vector 
r(r = r/r is a unit vector along r). The pair distribution function #(1,2) of the isotropic 
fluid is of particular interest as it is the lowest order microscopic quantity which contains 
informations about the translational and the orientational structures of the system and also 
has direct contact with intermolecular (as well as with intramolecular) interactions. For an 
ordered phase, on the other hand, most of the structural informations are contained in p(x) 
(see Sec. III). 

The values of the pair correlation functions as a function of intermolecular separation and 
orientations at a given temperature and pressure are found either by computer simulations 
or by solving the Ornstein-Zernike equation 



/i(l,2)-c(l,2)= 7(1,2) 

= p / |c(l,3)[ 7 (2,3)+c(2,3)]rf3 



(2.5) 
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where d3 = dr 3 dfls and h(l,2) = (7(1,2) — 1 and c(l,2) are, respectively, the total and 
direct pair correlation functions, using a suitable closure relation. Most commonly used 
closer relations are the Percus-Yevick (PY) and the hyper netted chain (HNC) relations. 
Approximations are introduced through these closure relations. The PY and HNC integral 
equation theories are given by the OZ equation coupled with the closure relation || 



0^(1,2) = /(l,2)[l + 7(1, 2)] (2.6) 

and 

C HNC (1, 2) = h{l, 2) - ln[l + h{l, 2)] - (3u(l, 2) (2.7) 

respectively. Here /(l, 2) = exp[— (3u(l, 2)] — 1 and (3 = (fe^T) -1 . 

Both the PY and HNC integral theories have been used to find the pair-correlations 
functions of model fluids of non-spherical molecules |], |10j- It is found that while the 
PY theory underestimates the correlations, particularly the angular correlation while the 
HNC theory overestimates them. In case of hard-core fluids we proposed a 'mixed' integral 
equation which interpolates between the HNC and PY theories and is thermodynamically 



consistent [IT]. Such an approach is needed for the soft-core potential the one considered in 
this paper also. We, however, defer this approach for the future and confine ourselves here 
to solve the PY equation to get the pair correlation functions for the GB(n-6) potential. 

The angle dependent function A(r 12 , f2i, f2 2 ) (where A may be pair correlation function 
or pair potential) is expanded in a basis set of rotational invariants || in space fixed (SF) 
frame according to the equation 



A(r ia ,n 1 ,n a ) = X; E ^^ 2 K^2)C g (Zii 2 i;m 1 m 2 m)^ imi (fi 1 )^ 2m2 (fi 2 )y;^(fi) (2.8) 

where C g{}\l 2 l\m\m 2 yn) are the Clebsch-Gordon coefficients. 

For fully axially symmetric particles it is also possible to expand the function in products 
of spherical harmonics in body fixed (BF) frame according to the equation. 



A(r 12 , fii, fi 2 ) = J2 A h i 2 m(n2)Y hm (^i)Y h! n(n 2 ) (2.9) 

l\l'2.m 
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where m = — m. Numerically it is easier to calculate the BF harmonic coefficients than 
the SF harmonic coefficients. The two harmonic coefficients are related through a linear 
transformation, 



In any numerical calculation we can handle only a finite number of the spherical harmonic 
coefficients for each orientation-dependent function. The accuracy of the results depends on 
this number. As the anisotropy in the shape of molecules (or in interactions) and the value 
of fluid density pj increases more harmonics are needed to get proper convergence. We have 
found that the series get converged if we truncate the series at the value of I indices equal 
to 6 for molecules with x < 3 11. Though it is desirable to include higher order harmonics 
i.e. for I > 6 but it will increase computational time many fold. Our interest is to use the 
data of the harmonics of pair correlation functions for freezing transitions where only low 
order harmonics are generally involved (see Sec. Ill below). The only effect the higher-order 
harmonics appear to have on these low-order harmonics is to modify the finer structure of 
the harmonics at small values of r whose contributions to the structural parameters (to be 
define below) are negligible. 

Using the numerical procedure outlined elsewhere ||, we have solved the PY equation 
for the GB(n-6) fluid having n values 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28 and 30 for 
x = 3.0 and well depth ratio k' — 5 at reduced temperatures, T* = ^ = 0.65, 0.80, 0.95 
and 1.25 for a wide range of densities. The other two parameters p and v are taken to be 
2 and 1, respectively. The solutions could be found only upto certain density p' the value 
of which depend upon the temperature and the value of n. The value of p' is often close to 
the isotropic-nematic transition. Because of this one faces problems in locating other less 
symmetric phases of the system using the theory to be discussed in Sec. III. 

In Fig. 2 we compare the values of g(r) = 1 + h ° ^ in BF frame at T* = 0.80 and density 
r](= |pj(TQa;o) = 0.25 for four sets of (n-6) combinations. It is seen from this figure that the 
first peak becomes sharper and attains its maximum value at smaller value of r*(— —) as 



the hardness of the core increases. The cause of this becomes clear if we look at Figs. 3 and 
4 which depict v (r) = — T* ln[(e~ /3 "^' ni ' n2 ^)n 1 ,Q 2 ] as a function of interparticle separation at 
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or 




in 



(2.10) 



T* = 0.8 andl.25, respectively. v(r) may be regarded as an averaged pair potential and, 
therefore, helps us in understanding the features of g(r). v(r) seems to have two minimum; 
one at r* rs 1.25 and other at r* rs 2.25. The first minimum becomes deeper at higher n and 
at lower temperature and almost vanishes at lower n and higher temperatures. The second 
minimum dependence on n (as well as on temperature ) is weak. One may also note the 
shift to lower values of r* of first minimum as n is increased. 

Since the PY theory is known to be reasonably accurate for systems interacting via pair 



potential which has hard repulsive core (n — > oo) and weak attraction [11], the values of the 
pair correlation functions reported here are expected to be more accurate for higher values 
of n and lower values of T* compared to values corresponding to lower n and higher T*. 

In Figs. 5-6 we compare the two other projections of PCF in BF- frame at the same state 
conditions and observe similar behavior. In Fig. 7 we compare the value of g(r) at r\ — 0.5 for 
GB(10-6) model at four different temperatures. Here we see that the first peak gets sharper 
as the temperature decreases. Such behavior is also seen (see Fig. 2) when n is increased at 
the same temperature. This is due to increasing tendency of the molecules to form parallel 
configurations. 

As has already been mentioned, the PY theory underestimates the molecular correlations. 
This can be seen from the pressure calculated using the values of the direct pair correlation 
function through the compressibility relation which is found to be lower than the simulated 
value as shown in Fig. 8. 

For a system consisting of axially symmetric non-dipolar molecules the static Kerr con- 
stant K is given by fl2], [13 



K = (5k 



I _ ~ 22 



5 

where C22 is structural parameter defined as in Eq.(3.22) and k is a constant dependent only 
upon single particle properties. The divergence of K may signal the absolute stability limit 



of the isotropic phase relative to orientationally ordered phase 1 12| . Thus the isotropic phase 
becomes orientationally unstable when the inverse Kerr constant K^ 1 — > 0. It is, however, 
important to emphasize that the condition K^ 1 — > does not determine the thermodynamic 
phase transition, but rather a point on the spinodal line. This means that the density at 
which K^ 1 = establishes a stability limit in the sense that at higher densities the isotropic 
phase cannot exist even as a metastable state. 

8 



The reduced Kerr constants (3AK 1 as a function of r\ for the various n values are plotted 
in Figs. 9 & 10 at T* = 0.8 and 1.25. 



III. THEORY FOR FREEZING 

The structural informations of fluids at the pair correlation functions level obtained above 
can be used to obtain information about their freezing. At the freezing point the spatial 
and orientational configurations of molecules undergo a modification. Often abrupt change 
in the symmetries of the system takes place on the freezing. In contrast to the isotropic 
fluid, the molecular configurations of most ordered phases are adequately described by the 
single particle density (singlet) distribution p(x).p(x) provides us with a convenient quantity 
to specify an arbitrary state of a system. One may consider a variational thermodynamic 
potential as a functional of p(x). The equilibrium state of the system at given T and P is 
described by the density p(T, P, x) corresponding to the minimum of the thermodynamic 
potential with respect to p(x). This forms the basis of the density functional theory. 

In this article we investigate the freezing of the GB(n-6) fluid into the nematic and the 
smectic A (Sm A) phases using density functional theory (DFT). In the nematic phase the 
full translational symmetry of the isotropic fluid phase (denoted as R 3 ) is maintained but 
the rotational symmetry 0(3) or SO (3) (depending upon the presence or absence of the 
centre of symmetry) is broken. In the simplest form of the axially symmetric molecules the 
group 0(3) (or SO (3)) is replaced by one of the uniaxial symmetry (or D^). The 

phase possessing the R 3 A (denoting the semi direct product of the translational group 
R 3 and the rotational group -Dooh) symmetry is known as uniaxial nematic phase [I], [14 . 



The smectic liquid crystals, in general, have a stratified structure with th long axes of 
molecules parallel to each other in layers. This situation corresponds to partial breakdown 
of translational invariance in addition to breaking of the orientational invariance. Since a 
variety of molecular arrangement are possible within each layer, a number of smectic phases 
are possible [|IJ. The simplest among them is the Sm A phase. In it the centre of mass of 
molecules in a layer are distributed as in a two-dimensional fluid but the molecular axes are 
on the average along a direction normal to the smectic layer (i.e. the director n is normal to 
the smectic layer). The symmetry of the Sm A phase is D OQ i l /\(R 2 x Z) where R 2 corresponds 
to a two-dimensional liquid structure and Z for a one-dimensional periodic structure. 
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The order parameters which characterize the ordered structures can be found from the 
singlet distribution p(x). For this we express it in the Fourier series and the Wigner rotation 
matrices. Thus 



p(x) = p(r, ft) = Po ]T £ Qimn(G q ) eMiG q .v)D l mn (n) (3.1) 

q Ivan 



where the expansion coefficients 

21 + 1 



Qimn(G q ) = Jdv J cMp(r, n) exp(-iG.r)£>£ n (n) (3.2) 

are the order parameters, G q the reciprocal lattice vectors, po the mean number density 
and D l mn (Q) the generalized spherical harmonics or Wigner rotation matrices |T3 . 

Since we are interested in uniaxial systems of cylindrically symmetric molecules, m = 
n = in Eqs.(3.1) & (3.2). This leads to 

p{r,n) =poEE^ ex p( iG ?- r ) p «( cos0 ) ( 3 - 3 ) 

i q 

and 

21 + 1 r r 

Qiq = J dv J dn P (r, n) expHG.r)P,(cos#) (3.4) 

where Pi(cos8) is the Legendre polynomial of degree I and 9 is the angle between the 
cylindrical axis of a molecule and the director. 

Since in the nematic phase the centres of mass of molecules are distributed as randomly as 
in the isotropic fluid but the molecular axes are aligned along a particular direction defined 
by the director n (a unit vector) we have G q = and 

Q m = ((21 + l)P,(cos(0))) = (21 + 1)P, (3.5) 

where angular bracket indicates the ensemble average. It is often enough to use two orien- 
tational order parameters P2 and P4 to locate the isotropic-nematic transition as in almost 
all known cases the transition is weak first-order transition [|TJ. 

To characterize the Sm A phase we need three different class of order parameters; (i) 
orientational, (ii) positional, (iii) mixed. These parameters are found from Eq.(3.2). For the 
orientational order we take P% and P4 as in case of the nematic phase. For the positional 
order along the z-axis we choose one order parameter corresponding G z = d being the 
layer spacing. Thus 
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2ttz 

/i = Qoo(Gy = (cos(— )) (3.6) 

The coupling between the positional and orientational ordering is described by the (mixed) 
order parameter r defined as 

1 2lTZ 

r = -Q 20 (G Z ) = (cos(— )P 2 (cos£)> (3.7) 
o a 



We therefore choose four order parameters to describe the ordering in a Sm A phase 
and two for the nematic ordering. Another way of writing the trial singlet distribution 
corresponding to the ordered phases of our interest is 

p(r, ft) = A p exp[-a{z - df -a\z- d) 2 P 2 (cos#) + A 2 P 2 (cosfl) + A 4 P 4 (cos#)] (3.8) 

where Aq is a normalization constant, a and a 1 are associated with the formation of layer 
in the Sm A phase and A 2 and A4 with orientational ordering. When a and a 1 are zero 
but A 2 and A 4 are non zero the phase is nematic. In case of the isotropic fluid all the four 
parameters a, a 1 , A 2 and A 4 are zero. If all the four parameters are non zero the phase is Sm 
A. The four order parameters defined above can be found taking the expression of p(r, ft) 
given by Eq.(3.8). Thus 

fj, — — / az cos{— — ) / dxexp{b) (3.9) 



d Jo d Jo 

Ao 
d 

Ao f d 
d 



rd 2,712 1'^ 

/ dzcos{— ) / dxexp(S)P 2 (x) (3.10) 
Jo a Jo 

f d dz f 1 dxexp(S)P 2 (x) (3.11) 
Jo Jo 



p 4 = _2 / dz I dxexp{S)P 4 {x) (3.12) 
d Jo Jo 

where S = -a(z - df - a x (z - d) 2 P 2 (x) + A 2 P 2 (a;) + A 4 P 4 (a;) 
A. Density Functional Approach 

In the usual density functional theory aapproach one uses the grand thermodynamic 
potential to locate the transition. The grand thermodynamic potential is defined as 
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-W = (3 A- (3p J dxp(x) 



(3.13) 



where A is the Helmholtz free energy, \x the chemical potential and p(x) is a singlet distri- 
bution function. Eq.(3.1) can be written as 



AW = W -W f = AW X + AW 2 
where Wf is the grand thermodynamic potential of the isotropic fluid, and [14 



(3.14) 



AW 

N 
and 
AW 2 



— f drdfl | p(r, fi) In 



PfV 



Pf 



Ap(r,n) 



(3.15) 



— / dT 12 dQ 1 dn 2 Ap(T 1 , ni)c(r la , ft 1: fl 2 )Ap(r 2 , Q 2 ) (3.16) 
lot J 



N 2p f 

Here Ap(x) = p(x) — p/, where pf is the density of the coexisting liquid. The ordered phase 
density is found by minimizing AW with respect to arbitrary variations in the ordered phase 
density subject to the constraint which corresponds to some specific features of the ordered 
phase. Thus, 



In 



Pf 



A L + y dv 2 dn 2 c(v 12 ,n 1 ,n 2 ;p f )Ap(r 2 ,n 2 ) 



(3.17) 



where Al is Lagrange multiplier which appears in the equation because of constraint imposed 
on the minimization. 

One attempts to find solution of p(x) of Eq.(3.17) which have symmetry of the ordered 
phase. These solutions, inserted in Eq.(3.14) give the grand thermodynamic potential dif- 
ference between the ordered and liquid phases. The phase with the lowest grand potential is 
taken as the stable phase. Phase coexistence occurs at the value of pf which makes ——jf- = 
for the ordered and liquid phases. Substituting Eq.(3.1), into Eq.(3.17) and Eq.(3.14) and 
integrating results in, respectively 



+ Wfl = VI dnd^e-^McosOj exp[A L + £ E ^e^C?^)] 

(3.18) 
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and 



AW 



where 



-Ap* + Ap*C ° )0 

I ^ \ V ■* QhqQh'q A,q 

+ 2fe^(2L + l)(2L' + l) L ' L ' 



(3.19) 



Qo.o = Ap* (3.20) 

c£ o (0i) = {2i + i) P , J dr 12 dn 2C {r 12 ,n 1 ,n 2 ) (3.21) 

e* G « ri2 P z (cos0 2 ) 

C* , = (2/ + 1)(2Z' + l)p f J dv 12 dSl 1 dSl 2 (3.22) 
e iG ' ri2 c(r 12 , n ls n 2 )P,(cos0 1 )P,/(cos0 a ) 

are the structural parameters related to the Fourier transformed direct correlation function 
of the fluid phase. Eq.(3.18) is the expression for the order parameters. This version of the 
density functional theory is known as the second order density functional (SODFT) because 
it considers only the pair correlation functions and neglects the higher order correlations 
which might be present in the system at the transition point. 



B. Modified Weighted-Density Approximation (MWDA) 

In another version of the density functional approach in which higher order correlations 
are included and known as Modified Weighted Density Approximation PB[, one uses the 



Helmholtz free energy to locate the transition. For the Helmholtz free energy we write 



A[p{r, □)] = A id [p{r, SI)] + A ex [p(v, SI)} (3.23) 

where both terms in Eq.(3.23) are unique functionals of the one-particle density p(r,Si). 
The first term in the right hand side of Eq.(3.23) is a non uniform ideal gas contribution of 
the form 



A ld [p(r, SI)] = /T 1 j drdSlp(v, Sl){\n[p(r, Sl)\ 3 ] - 1} (3.24) 
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where A is the thermal de Broglie wavelength. The second term in the right hand side of 
Eq.(3.23) is the excess Helmholtz free energy of the non uniform system. 

In the modified weighted-density approximation the excess free energy of a uniform sys- 
tem, but evaluated of a weighted density p WE\ 



A MWDA [p] = {3 25) 

where N is the number of particles in the system 4>o(p) is the excess free energy per particle 
of a uniform system at density p. The weighted density p is constructed from the actual 
inhomogeneous one-particle density p(x) and is defined by 



P N 



^ J v dxp(x) J v dx'p(x')d;(x - x'; p) (3.26) 

introducing thereby the weighted function u)(x — x';p). It is an essential ingredient of 
the MWDA that the weighted function u which is used to determine the weighted density, 
depends itself on the sought function p; thus Eq.(3.26) has to viewed as a self-consistency 
condition for the determination of the weighted density. To ensure that the approximation 
in the determination of p becomes exact in the uniform limit, the weighted function has to 
be normalized, i.e., 

J dx£(x - x; p) = 1 (3.27) 

for any p. The function uj can be then uniquely specified by requiring that the approximate 
functional A^ c WDA [p] is exact upto second order in the functional expansion, namely 



C(x — x'; p ) = —3 lim 

P^PO 



5 2 A MWDA [p] 



^3.28) 



5p(x)5p(x') 

The conditions [Eqs. (3.25-3.28)] result in a particularly simple expression for uj, namely 



u;(x-x';p) 



^C(x-x';p) + lp0o(p) 



(3.29) 



20^(p) 

where V is the volume of the sample, <po{p) is the excess free energy per particle of an 
isotropic fluid of density p and primes on 0o(p) indicate derivatives with respect to density. 
Using expansion (Eq. 3.1) and (Eq.2.8), respectively for p(x) and C(x — x'; p) we find for 
the ordered phase 
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9 ~ P ° 5 5 ? (2L! + 1K2L 2 + 1) 



2pM(p) 



PoP 



0o (p) 
20 O (P) 



(3.30) 



Having computed p, the next step in freezing analysis is to substitute p into Eq.(3.25) to 
compute A^, WDA . In terms of structural parameter, the excess free energy per particle of a 
uniform system at a density p is given as 



PMp) 



Jo p"Z Jo 



(3.31) 



The ideal gas part is calculated using the ansatz for p(r, ft) given by Eq.(3.8). Thus 



PA id [p ] = Po J dr^J2^==^ Gq - r [{HAoPo\ 3 ) - 1}6lo ~ a(z - x ) 2 ^ 

+ + — g a (z - x ) 5lo] (3-32) 

To determine the transition parameters, we first compute the effective density p from 
Eq.(3.30) and minimizing the free energy from Eqs.(3.23, 3.31 and 3.32) with respect to 
p , «, A 2 , A 4 and a 1 . In order to determine the transition density of the coexisting isotropic 
(pf) and anisotropic (p ) phases it is necessary to equate the pressure and chemical potentials 
(Maxwell construction) of the two phases. 



IV. RESULTS AND DISCUSSION 



We have used both versions of the density functional methods described above to locate 
the freezing transitions and calculate the values of the freezing parameters. The structural 
parameters defined by Eq.(3.22) which appear in the density functional theory as the input 
data are obtained from the harmonics of the direct pair correlation functions evaluated 
using the PY integral equation theory (given in Sec. II). Using these values of the structural 
parameters and the four order parameters P2,Pa,p andr we have solved Eqs.(3.18-3.19) of 
the SODFT and Eqs. (3.23-3.32) of the MWDA for the GB(n-6) fluid with 8 < n < 30 for 
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temperatures lying between 0.65 to 1.25. All our results correspond to /i = 2, v — 1, xo = 
3 and k' — 5. 

Our results show that for none of the cases studied here Sm A phase gets stabilized. In 
the low temperature region for a given n it, however, appeared as a metastable state having 
free energy lower than that of the isotropic phase but higher than the nematic (see Table I). 
Since we have not included Sm B and crystalline phases in our investigation for the reason 
already given, we found only the isotropic-nematic transition. 

For each n we found a lower cut-off of the temperature for the existence of the nematic 
phase. The nematic phase was not found to exist below this temperature. The lower cut-off 
temperature for the nematic phase is found to increase with n. For example, where for 
n = 8 and 10 we found the nematic phase to exist at T* = 0.65 but not for n > 12. The 
computer simulation results of Miguel et al |J show that for n — 12 the cut-off temperature 
is slightly above T* = 0.8. Our results, however, show that the nematic phase exists at 
T* = 0.8. This may be due to error in the structural parameters values found from the PY 
theory. 

Both versions of the density functional theory give similar results for the transition den- 
sity p*f but give the values of the order parameters including the change in density at the 
transition which are different from each other. More surprising is the way the values of 
the order parameters Pi and P4 vary with temperature and with n (see Tables II- V) found 
from the two version of the theory. While the SODFT predicts that Pi and P4 decrease as 
the transition temperature is increased, the MWDA predicts them to increase. The com- 
puter simulation results || |6j do not give any clear indication as how these parameters vary 
with transition temperature. Similar difference in the variation of the values of the order 
parameters with n is also found. 

In Table II- V we give the values of the transition parameters found from the two theories. 
We also give the results found from the computer simulations at T* = 0.95 and 1.25 
for n = 12. There is very good agreement between these results at T* = 0.95. The 
transition density found from the theories are identical though somewhat higher than the 
value found from the simulation. The value of Ap* found from these methods are also in 
good agreement, though MWDA predicts the value of Ap* which is lower than the SODFT 
as well as simulation value. Pressure and chemical potentials are in good agreement. But 
there is difference in the value of Pi and P4. At T* = 1.25 both theories predict the transition 
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density which is high compared to the MD value. One of the possible reasons for this is, as 
pointed out in Sec. II, the inaccuracy in the values of c-harmonics at higher temperature. 
The PY theory is known to underestimate the angular correlations and this defect of the PY 
theory becomes more pronounced as temperature is increased for a given n. This may be 
the reason why the theory predicts the transition at higher density than the MD value. As a 
consequence of this the transition pressure and the chemical potential are also substantially 
higher than the MD values. This comparison at T* = 0.95 and 1.25 show that while the 
DFT is good to predict the freezing parameters, the PY values of structural parameters at 
higher temperature are lower than the actual values. 

We hope to combine the PY and HNC theories to generate accurate values of the har- 
monics of the pair correlation functions and with these values to compute the full phase 
diagram. 
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TABLE I: Values of order parameters and energy of smectic A and nematic phases at T* = 0.8 
for GB(12-6) potential. While nematic is a stable phase, smectic A is metastable as its energy is 
higher than the nematic 



* 

pf 


Phase 




Pz 




h 


Pa 


T2z 




Ap* 


AW 


0.293 


Sm-A 




0.674 




0.891 


0.719 


0.687 




0.098 


-0.023 




Nematic 




0.000 




0.843 


0.562 


0.000 




0.058 


-0.093 


0.306 


Sm-A 




0.647 




0.888 


0.741 


0.677 




0.088 


-0.157 




Nematic 




0.000 




0.903 


0.654 


0.000 




0.063 


-0.243 


0.312 


Sm-A 




0.629 




0.887 


0.750 


0.667 




0.083 


-0.234 




Nematic 




0.000 




0.923 


0.691 


0.000 




0.065 


-0.333 


TABLE II: Isotropic- Nematic transition parameters 


for GB(n- 


•6) fluid at T* = 


0.65. The reduced 


units are P* 


= Pa$/e ,p* = 


: /V e o, 


and 


* 3 
P = P°0 












Potential Model 


Theory 


Pf 


# 
Pn 


Ap* 




Pa 


p* 


* 

P 


(8,6) 




DFT 




0.428 


0.431 


0.006 


0.69 


0.40 


9.48 


26.79 






MWDA 


0.412 


0.416 


0.009 


0.56 


0.27 


7.77 


22.71 


(10, 6) 




DFT 




0.29 


0.301 


0.038 


0.72 


0.41 


1.30 


3.96 






MWDA 


0.286 


0.29 


0.013 


0.36 


0.12 


1.22 


3.69 
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TABLE III: Same as in Table II but at T* = 0.80 



Potential Model 


Theory 


* 

Pf 


* 

Pn 


Ap* 


P2 




p* 


* 

P 


(10, 6) 


DFT 


0.341 


0.346 


0.015 


0.68 


0.37 


3.81 


12.44 




MWDA 


0.339 


0.341 


0.007 


0.40 


0.15 


3.71 


12.16 


(12, 6) 


DFT 


0.282 


0.295 


0.046 


0.74 


0.43 


1.38 


4.11 




MWDA 


0.277 


0.281 


0.015 


0.36 


0.12 


1.27 


3.69 


(Id fi") 


DFT 


0.239 


977 


1 58 


99 


69 


58 


n 65 




MWDA 


0.237 


0.241 


0.019 


0.27 


0.08 


0.55 


0.54 




TABLE IV: Same 


as in Table 


II but a 


t T* = 0.95 








Potential Model 


Theory 


* 

Pf 


* 

Pn 


Ap* 




Pa 


p* 


* 

P 


(10, 6) 


DFT 


0.381 


0.385 


0.009 


0.68 


0.38 


8.26 


25.63 




MWDA 


0.379 


0.382 


0.007 


0.46 


0.20 


8.04 


25.04 


(12, 6) 


MD 


0.308 


0.314 


0.019 


0.50 




3.50 


12.70 




DFT 


0.322 


0.328 


0.02 


0.67 


0.37 


3.40 


11.28 




MWDA 


0.322 


0.325 


0.008 


0.37 


0.13 


3.40 


11.28 


(14, 6) 


DFT 


0.287 


0.299 


0.042 


0.74 


0.43 


1.82 


5.66 




MWDA 


0.283 


0.288 


0.017 


0.37 


0.12 


1.69 


5.21 


(16, 6) 


DFT 


0.261 


0.283 


0.085 


0.82 


0.51 


1.06 


2.59 




MWDA 


0.245 


0.251 


0.027 


0.36 


0.12 


0.82 


1.64 
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TABLE V: Same as in Table II but at T* = 1.25 



Potential Model 


Theory 


* 

Pf 


* 

Pn 


Ap* 


h 


Pa 


p* 


* 

A* 


(10, 6) a 


DFT 


0.454 


0.456 


0.005 


0.72 


0.44 


26.93 


72.59 




MWDA 


0.435 


0.437 


0.005 


0.62 


0.30 


21.26 


59.84 


(12, 6) 


MD[5] 


0.323 


0.331 


0.025 


0.50 




5.70 


20.90 




DFT 


0.378 


0.382 


0.009 


0.68 


0.38 


10.90 


34.27 




MWDA 


0.375 


0.378 


0.007 


0.47 


0.21 


10.42 


32.99 


(14, 6) 


DFT 


0.344 


0.349 


0.014 


0.68 


0.38 


6.57 


21.88 




MWDA 


0.343 


0.346 


0.009 


0.44 


0.20 


6.52 


21.71 


(18, 6) 


DFT 


0.306 


0.315 


0.028 


0.72 


0.41 


3.43 


11.52 




MWDA 


0.303 


0.307 


0.014 


0.41 


0.18 


3.25 


10.95 


(24, 6) 


DFT 


0.273 


0.291 


0.065 


0.79 


0.49 


1.81 


5.35 




MWDA 


0.267 


0.274 


0.026 


0.37 


0.13 


1.65 


4.78 


(30, 6) 


DFT 


0.249 


0.283 


0.137 


0.90 


0.61 


1.11 


2.38 




MWDA 


0.242 


0.248 


0.027 


0.34 


0.11 


0.99 


1.90 



°The results have been found by extrapolating the data of the structural parameters to high densities. The 
value of the transition parameters may, therefore, not be as accurate as for the other cases. 
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